First I need to load up the packages I’ll need
library(sf)
library(ggplot2) #development version!
## devtools::install_github("tidyverse/ggplot2")
library(tidyverse)
library(readr)
## Not sure about this bit
#library("tidyverse",lib.loc="/Library/Frameworks/R.framework/Versions/3.4/Resources/library")
library(cowplot)
library(sp)
library(gridExtra)
library(dplyr)
Now I import my data. I filter for the Arran postcodes, (since Arran all begins ‘KA27’).
## Finding the Arran coordinates
allcoordinates <- read.csv("alldata/ukpostcodes.csv")
arrancoordinates <- filter(allcoordinates,substr(postcode,1,4)=="KA27")
Now I plot these coordinates.
## Plotting the Arran coordinates
ggplot(data = arrancoordinates) +
geom_point(mapping = aes(x = longitude, y = latitude), shape=20) +
ggtitle("Arran Postcodes") +
labs(title = "Arran Postcodes", x = "Longitude", y = "Latitude") +
theme(plot.title = element_text(hjust = 0.5)) +
theme_grey() +
coord_map()
Now I create some plots. #Arran Borders
pcs <- read_sf("alldata/Scotland_pcs_2011")
#Print Post codes lists
arransubsect <- filter(pcs,substr(label,1,4)=="KA27")
arransubsect %>%
ggplot() +
geom_sf() +
theme(axis.text.x=element_text(angle=45, hjust = 1)) +
theme_grey()
After a little editing I can overlay the two.
simple.sf <- st_as_sf(arrancoordinates, coords=c('longitude','latitude'))
st_crs(simple.sf) <- 4326
simple.sf %>%
ggplot() +
geom_sf(shape=20) +
theme_grey()
arransubsect %>%
ggplot() +
geom_sf() +
theme(axis.text.x=element_text(angle=45, hjust = 1)) +
theme_grey() +
geom_sf(data=simple.sf, shape=20)
Now I load the SIMD data, containing the geometries (shapefiles) and SIMD data (percentiles, etc)
#Import SIMD data from http://www.gov.scot/Topics/Statistics/SIMD
#https://data.gov.uk/dataset/scottish-index-of-multiple-deprivation-simd-2012
#https://data.gov.uk/dataset/scottish-index-of-multiple-deprivation-simd-2012/resource/d6fa8924-83da-4e80-a560-4ef0477f230b
DZBoundaries2016 <- read_sf("./alldata/SG_SIMD_2016")
DZBoundaries2012 <- read_sf("./alldata/SG_SIMD_2012")
DZBoundaries2009 <- read_sf("./alldata/SG_SIMD_2009")
DZBoundaries2006 <- read_sf("./alldata/SG_SIMD_2006")
DZBoundaries2004 <- read_sf("./alldata/SG_SIMD_2004")
I have to choose the right columns manually in order to select the Arran data.
#Selecting Arran data from Scotland (2016)
#Find postcode look-up from below file for KA27 postcodes. Find unique DZ. Find row positions.
#SIMD2016 <-read.csv("./alldata/00505244.csv")
#Selecting ArranDZ2016
Arrandz2016 <- c(4672,4666,4669,4671,4667,4668,4670)
arran2016 <- DZBoundaries2016[Arrandz2016,]
#Find postcode look-up, KA27 postcodes. Find unique DZ. Find row positions.
#Selecting ArranDZ2012
Arrandz2012 <- c(4409,4372,4353,4352,4351,4350,4349)
#2012
arran2012 <- DZBoundaries2012[Arrandz2012,]
#2009
arran2009 <- DZBoundaries2009[Arrandz2012,]
#2006
arran2006 <- DZBoundaries2006[Arrandz2012,]
#2004
arran2004 <- DZBoundaries2004[Arrandz2012,]
Now I want to plot all the data, first I combine it all into one table. First I subselect the data I want from the appropriate columns.
arran20162 <- arran2016 %>%
select(DataZone, geometry, Percentile) %>%
mutate(year="2016")
arran20122 <- arran2012 %>%
select(DataZone, geometry, Percentile) %>%
mutate(year="2012")
arran20092 <- arran2009 %>%
select(DataZone, geometry, Percentile) %>%
mutate(year="2009")
arran20062 <- arran2006 %>%
select(DataZone, geometry, Percentile) %>%
mutate(year="2006")
arran20042 <- arran2004 %>%
select(DataZone, geometry, Percentile) %>%
mutate(year="2004")
arransimd <- rbind(arran20162,arran20122,arran20092,arran20062,arran20042)
The reason I’ve downloaded all the datazones shapefiles individually is because they change between 2016 and 2012. See the small differences.
arran1612 <- rbind(arran20162,arran20122)
arran1612 %>%
ggplot() +
geom_sf(aes(fill = DataZone)) +
facet_wrap('year') +
theme_grey() +
theme(legend.position="none") +
theme(axis.text.x=element_text(angle=45, hjust = 1))
Now I plot the percentiles.
arransimd %>%
ggplot() +
geom_sf(aes(fill = Percentile)) +
facet_wrap('year') +
theme_grey() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(), axis.text.y=element_blank(),
axis.ticks.y=element_blank())
arransimd %>%
ggplot() +
geom_sf(aes(fill = Percentile)) +
facet_wrap('year', nrow = 1) +
theme_grey() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(), axis.text.y=element_blank(),
axis.ticks.y=element_blank())
There we are. Not the SIMD health percentiles of Arran zones throughout SIMD history. And I’ve learned a little bit about graphics in R.
If I wanted to I could show the zones individually.. First I find the unique zones. (There are 14. 7 Zones 2016, 7 Zones pre-2016)
datazones <- unique(arransimd$DataZone)
datazones
[1] "S01011177" "S01011171" "S01011174" "S01011176" "S01011172" "S01011173" "S01011175" "S01004409" "S01004372" "S01004353"
[11] "S01004352" "S01004351" "S01004350" "S01004349"
I’ll have to find out a simpler way to do this but.. #Pre-2016 Individual Zones
function0.5 <- function(argument)
{
filter(arransimd, DataZone==argument)
}
filter(arransimd, DataZone=="S01004353")
Simple feature collection with 4 features and 3 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 186117 ymin: 618919 xmax: 206903.5 ymax: 652363.5
epsg (SRID): NA
proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
DataZone Percentile year geometry
1 S01004353 78 2012 MULTIPOLYGON(((200930 63676...
2 S01004353 68 2009 MULTIPOLYGON(((200930 63676...
3 S01004353 71 2006 MULTIPOLYGON(((200930 63676...
4 S01004353 65 2004 MULTIPOLYGON(((200930 63676...
function0.5("S01004353")
Simple feature collection with 4 features and 3 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 186117 ymin: 618919 xmax: 206903.5 ymax: 652363.5
epsg (SRID): NA
proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
DataZone Percentile year geometry
1 S01004353 78 2012 MULTIPOLYGON(((200930 63676...
2 S01004353 68 2009 MULTIPOLYGON(((200930 63676...
3 S01004353 71 2006 MULTIPOLYGON(((200930 63676...
4 S01004353 65 2004 MULTIPOLYGON(((200930 63676...
pre2016list2 <- list("S01004409", "S01004372", "S01004353", "S01004352", "S01004351", "S01004350", "S01004349")
#Create a new way of making character list
pre2016list <- lapply(pre2016list2, function0.5)
function1 <- function(argument)
{
argument %>%
ggplot() +
geom_sf(aes(fill = Percentile)) +
facet_wrap('year') +
theme_grey() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
}
function2 <- function(argument)
{
arransubsect %>%
ggplot() +
geom_sf() +
theme_grey() +
theme(axis.text.x=element_text(angle=45, hjust = 1)) +
geom_sf(data= argument, aes(fill = DataZone))
}
function3 <- function(argument)
{
argument %>%
ggplot() +
geom_sf(data = arransubsect) +
geom_sf(aes(fill = Percentile)) +
facet_wrap('year') +
theme_grey() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
}
function4 <- function(argument)
{
plot_grid((function1(argument)), (function2(argument)), labels = c("A", "B"))
}
lapply(pre2016list, function4)
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
Now I apply the same functions to the post-2016 Zones
arransubsect <- filter(pcs,substr(label,1,4)=="KA27")
#post2016list <- list(S01011177, S01011171, S01011174, S01011176, S01011172, S01011173, S01011175)
post2016list2 <- list("S01011177", "S01011171", "S01011174", "S01011176", "S01011172", "S01011173", "S01011175")
post2016list <- lapply(post2016list2, function0.5)
lapply(post2016list, function4)
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
Pre-2016
arransimd2 <- dplyr::filter(arransimd, year < 2016)
arransimd2 %>%
ggplot(aes(x=year, y=Percentile)) +
geom_bar(stat="identity") +
facet_wrap('DataZone') +
theme_grey() +
labs(title = "Arran SIMD Datazones", x = "Year", y = "Percentile") +
theme(plot.title = element_text(hjust = 0.5))
Post-2016
arransimd3 <- dplyr::filter(arransimd, year >= 2016)
arransimd3 %>%
ggplot(aes(x=year, y=Percentile)) +
geom_bar(stat="identity") +
facet_wrap('DataZone') +
theme_grey() +
labs(title = "Arran SIMD Datazones", x = "Year", y = "Percentile") +
theme(plot.title = element_text(hjust = 0.5))
Now I want to overlay the postcodes for a particular shapefile, in this case by Datazone. To do this I’ve converted both the Arran coordinates and Arran (2016) shapefiles into Spatial Points/Polygons, converted them into a common CRS, and then compared them by using over().
exampleshapes <- sf:::as_Spatial(arran2016$geom)
examplepoints <- sf:::as_Spatial(simple.sf$geom)
examplepoints <- spTransform(examplepoints, CRS("+proj=longlat +datum=WGS84"))
exampleshapes <- spTransform(exampleshapes, CRS("+proj=longlat +datum=WGS84"))
namingdzpostcode <- over(exampleshapes, examplepoints, returnList = TRUE)
I can then take a member reference from the orginal postcode list, which gives me a selection of the rows in that DZ. For simplicity I’ve written this as a new function.
Unfortunately, I haven’t worked out how to coordinate the new ID with the original DZ names yet, so I have to select by using the appropriate ID for each DZ.
pre2016listID <- list(3,2,1,4,7,6,5)
post2016listID <- list(1,2,3,4,5,6,7)
function6 <- function(argument)
{
simple.sf[namingdzpostcode[[argument]],]
}
I can then use the above, and plot over the appropriate DZ shapefile. e.g #Projecting the coordinate selections #need to automate
function1(function0.5("S01004372")) +
geom_sf(data=function6(2), shape=20)
function2(function0.5("S01004372")) +
geom_sf(data=function6(2), shape=20)
function3(function0.5("S01004372")) +
geom_sf(data=function6(2), shape=20)
If I edit function 4 a little so that the geom_sf layer is a second argument then I can also use function 4.
function4.5 <- function(argument, argument2)
{
a <- function1(argument)
b <- function2(argument) +
geom_sf(data=function6(argument2), shape=20)
plot_grid(a, b, labels = c("A", "B"))
}
function4.5(function0.5("S01004372"),2)
I’ve also made another function to plot the DZ on it’s own with coordinates.
function5 <- function(argument, argument2)
{
argument %>%
ggplot() +
geom_sf() +
theme_grey() +
geom_sf(data=function6(argument2), shape=20) +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
}
function5(function0.5("S01004372"),2)
function1.5 <- function(argument)
{
argument %>%
ggplot() +
geom_sf(aes(fill = Percentile)) +
facet_wrap('year') +
theme_grey() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
theme(legend.position="bottom")
}
function2.5 <- function(argument)
{
arransubsect %>%
ggplot() +
geom_sf() +
theme_grey() +
theme(axis.text.x=element_text(angle=45, hjust = 1)) +
theme(legend.position="none") +
geom_sf(data= argument, aes(fill = DataZone))
}
function7.5 <- function(argument, argument2)
{
a <- function1.5(argument)
b <- function2.5(argument)
c <- function5(argument, argument2)
grid.arrange(a, b, c, nrow = 1)
}
function2.5.1 <- function(argument)
{
arransubsect %>%
ggplot() +
geom_sf() +
theme_grey() +
theme(axis.text.x=element_text(angle=45, hjust = 1)) +
theme(legend.position="bottom") +
geom_sf(data= argument, aes(fill = DataZone))
}
function7.5.1 <- function(argument, argument2)
{
a <- function1.5(argument)
b <- function2.5.1(argument)
c <- function5(argument, argument2)
grid.arrange(a, b, c, nrow = 1)
}
function7.5.1.2 <- function(argument)
{
function7.5.1(pre2016list[[argument]],pre2016listID[[argument]])
}
lapply(1:7, function7.5.1.2)
[[1]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[2]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[3]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[4]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[5]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[6]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[7]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
function7.5.1.1 <- function(argument)
{
function7.5.1(post2016list[[argument]],post2016listID[[argument]])
}
lapply(1:7, function7.5.1.1)
[[1]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[2]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[3]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[4]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[5]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[6]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
[[7]]
TableGrob (1 x 3) "arrange": 3 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]
plot text of percentile, at centre of shape file coordinates. overlay postcode labels.
overlay info onto leaflet then with labels.